/*
 * parmerge.c
 *
 *  Created on: Mar 12, 2019
 *      Author: sascha
 */

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#include "merge.h"

int rank(double x, double X[], int n) {
  int L = 0;
  int R = n;
  while (L < R) {
    int i = (L + R) / 2;
    if (X[i] < x) {
      L = i + 1;
    } else {
      R = i;
    }
  }
  return L;
}

int rank_non_strict(double x, double X[], int n) {
  int L = 0;
  int R = n;
  while (L < R) {
    int i = (L + R) / 2;
    if (X[i] <= x) {
      L = i + 1;
    } else {
      R = i;
    }
  }
  return L;
}

void merge(double a[], long n, double b[], long m, double c[]) {
    // replace this by a parallel merge algorithm
    //seq_merge1(a, n, b, m, c);

    int chunk_size;
    #pragma omp parallel
    {
        #pragma omp master
        chunk_size = n/omp_get_num_threads();
    }

    #pragma omp parallel for schedule(static, chunk_size)
    for (int i = 0; i < n; i++) {
      c[i + rank(a[i], b, m)] = a[i];
    }

    #pragma omp parallel for schedule(static, chunk_size)
    for (int j = 0; j < m; j++) {
      //c[j + rank_non_strict(b[j], a, n)] = b[j];
      c[j + rank(b[j], a, n)] = b[j];
    }

    //#pragma omp parallel for schedule(static, chunk_size)
    //for (int i = 0; i < n + m; i++) {
        //if (i < n)
	    //c[i + rank(a[i], b, m)] = a[i];
	//else if (i < n + m) {
	    //int j = i - n;
	    //c[j + rank(b[j], a, n)] = b[j];
	//}
    //}
}
